{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## FILL IN YOUR NAME AND UTORID HERE\n", "\n", "Python Assignment 2: Orbital motion \n", "==========================================================\n", "\n", "## Before you begin\n", "\n", "Assuming you have loaded this file into your Jupyter Notebooks workspace, make sure to press the \"play\" button at the top of the page in each box. This will render the Markdown into nicely-formatted text and execute all of the sections of Python code. Note that some sections will give error messages until you add the calculations you are asked for. \n", "\n", "## Introduction\n", "\n", "In class we have been discussing gravitation and the orbits of planets around the Sun. In this python notebook you will show how a planet's orbital speed changes over the course of its orbit, for circular orbits and for elliptical orbits. \n", "\n", "For this assignment, it will be useful to remind yourself of C&O sections 2.2 and 2.3. \n", "\n", "As you modify each of the following segments of code, remember to hit the \"play\" button at the top of the notebooks interface so that the cell containing the code \"runs\" (you can also hit Shift+Return for the same functionality). It should display the output of the piece of code just below that cell. You can run every cell in this notebook before you modify it, just to see what it produces. If you find any of this confusing or overwhelming, please refer back to the programming resources in the \"Guides\" section on Quercus. The \"Introduction to Programming\" guide contains links to some documents that should help you." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Here we will import the usual libraries\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "%matplotlib inline\n", "\n", "# Here we are adding a couple new libraries. These are extremely useful for keeping track of \n", "# units and constants! \n", "# Remember that when we import a library as 'x', we can then call all that library's functions by x.function()\n", "from astropy import units as u\n", "from astropy import constants as c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the astropy libraries to track constants and units ##\n", "Below is some simple code showing you how to use the astropy.units and astropy.constants libraries to better keep track of your calculations. In both coding and problem sets, if you get a result that seems off, check your units! \n", "\n", "For more information on astropy units, check out the [documentation](https://docs.astropy.org/en/stable/units/). " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$925 \\; \\mathrm{\\frac{km\\,s}{h}}$" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Play around with this to see how units work!\n", "# Simple calculation - constant motion in a straight line\n", "\n", "# Set a speed:\n", "v = 25. *u.km/u.hour\n", "# Set a time:\n", "t = 37 * u.s\n", "# Calculate distance traveled:\n", "d = v * t\n", "d" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$0.25694444 \\; \\mathrm{km}$" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Well that's not so useful! You can see that we have two units of time in the result. \n", "# Try this:\n", "d.to(u.km)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "My value for d without specifying units is 925.0\n", "My value for d in units of km is 0.2569444444444444\n", "My value for d in units of au is 1.7175675244717535e-09\n" ] } ], "source": [ "# You can also simply output the value of the result, but be sure to specify your desired unit or your answer\n", "# may be nonsensical\n", "print('My value for d without specifying units is {0}'.format(d.value))\n", "print('My value for d in units of km is {0}'.format(d.to(u.km).value))\n", "print('My value for d in units of au is {0}'.format(d.to(u.au).value))\n", "\n", "# Also this is way too many sig figs. Below I will include an example of formatted output. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up the problem ##\n", "\n", "There are a number of parameters that are needed to fully describe the motion of a planet around a star, or of two stars orbiting a common centre of mass. Avoiding the effects of orbital inclination and tilt, they are:\n", "- the semi-major axis $a$\n", "- the eccentricity $e$\n", "- the masses of the objects $m_1$ and $m_2$, or the mass of the star $M$ if $M >> m$\n", "- the period of the orbit $P$\n", "\n", "From these, we can calculate the distance $r$ of the planet from the star, and the orbital velocity $v$. \n", "\n", "Use the next block to set up your orbital parameters. To start, we will look at the orbit of an Earth-like planet in an Earth-like orbit around a star with the same mass as the Sun. " ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": true }, "outputs": [], "source": [ "## Question 1: Modify the code below to set up the problem as described above. Use the parameters for Earth's orbit. \n", "p = 427.*u.day # Use astropy units library to keep things nice\n", "a = 4.3*u.au # \n", "e = 0.4 # remember e has no units\n", "M = 1.6*u.Msun # Mass of the star\n", "m = 2.3*u.Mearth # Mass of the planet" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculating r and v\n", "\n", "When we assume that the orbit is circular, the distance from the planet to the star, and its orbital speed, remain constant. Here, we want to plot how the distance from the planet to the star, and the orbital speed, change over the course of the planet's orbit in the elliptical case. \n", "\n", "You will want to recall two equations:\n", "\n", "$r = \\frac{a(1-e^2)}{1 + e\\mathrm{cos}\\theta}$ \n", "\n", "$v = \\sqrt{GM (2/r - 1/a)}$\n", "\n", "You can see from C&O 2.3 that this solution for $v$ is for the motion of the centre of mass of a two body system, and the above relationship holds when $M >> m$. \n", "\n", "Refer to Figure 2.4 in C&O to remind yourself that $\\theta$ is the angle that describes the planet's position in its orbit, measured counterclockwise from the major axis of the ellipse. To plot $r$ and $v$, we need $\\theta$. " ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([3. , 3.33333333, 3.66666667, 4. , 4.33333333])" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## QUESTION 2: Create an array for theta in radians that will allow us to plot r and v over \n", "# *the entire orbit*. \n", "# We are using the np.linspace function, where:\n", "# array = np.linspace(start_value, stop_value, number of steps)\n", "# This will create an array of values between your start and stop value, evenly spaced with an interval between \n", "# values that depends on the number of steps you choose. \n", "# Replace these values the start and stop values you would need to plot a full orbit. \n", "# The number of steps to use is up to you; if your plots below are choppy you may want to go back and increase the number. \n", "\n", "theta_rad = np.linspace(3, 6, 10)\n", "\n", "# Print out the first few numbers:\n", "theta_rad[0:5]" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 9)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m9\u001b[0m\n\u001b[0;31m r =\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "## QUESTION 3: Write out code that will calculate r and v following the above equations, using your theta_rad array. \n", "# Remember that you can do calculations with arrays the same way you can do calculations with single numbers!\n", "# Use the variables a, e, M, etc. from above in your calculation, as we did in the simple\n", "# distance calculation above. \n", "# You divide with the '/' symbol, multiply with '*'. Use brackets to ensure the math is done in the correct order!\n", "# In python, power laws are calculated via x**2. \n", "# You can use np.cos(x) for cosine (assuming x is in radians), and np.sqrt(x) for square roots. \n", "# You can also do roots via power laws, like (x)**(1/2.), as an example. \n", "\n", "r = \n", "\n", "v = \n", "\n", "# Print out the first r and v value:\n", "# Note how we can adjust the decimal places being reported! \n", "# If you have used units, they will be included in your output. \n", "# If you haven't, this will produce an error! \n", "print('At r = {0:.4}, the velocity of the planet is {1:.4}'.format(r.to(u.au)[0], v.to(u.km/u.s)[0]))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# SCRATCH BOX\n", "# The box above will throw errors if your units aren't correct for r and v. \n", "# If you're having trouble making your units work, use this box as a scratch area\n", "# to test your calculations. You can always make a scratch box on your own in these \n", "# templates, just label it like this one, or delete it when you've solved your problem.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the results" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '')" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAD8CAYAAAA2edOFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEHZJREFUeJzt3N+LpXd9B/D3p7uGIlqizUbX3U03tXvRbSkYhhCwF9IY2awh60UvErAGvVgCDSit2LX5B/wBVUKDIVghwZQgqLjISoxpb2MziSaSrjFr0Gbd1axeqJCLsPjpxZwtk/Hszpk9Z2a+aV4vOJzzPM/3ec47M+GTd845c6q7AwDA9vqD7Q4AAIBSBgAwBKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYwM7tDnA5rrrqqt6/f/92xwC20JNPPvnL7t613TnmZX7B68+s8+s1Wcr279+f5eXl7Y4BbKGq+ul2Z1gE8wtef2adX96+BAAYgFIGADAApQwAYABKGQDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADGAhpayqDlXVc1V1qqqOTTleVXXP5PgzVXXdmuM7qup7VfXNReQB2AgzDBjB3KWsqnYkuTfJzUkOJrm9qg6uWXZzkgOT29EkX1hz/KNJTs6bBWCjzDBgFIt4pez6JKe6+4XufiXJw0mOrFlzJMmDveLxJFdW1e4kqaq9Sd6f5IsLyAKwUWYYMIRFlLI9SV5ctX16sm/WNZ9P8okkv1tAFoCNMsOAISyilNWUfT3Lmqq6JclL3f3kuk9SdbSqlqtq+dy5c5eTE2CaTZ9h5hcwi0WUstNJ9q3a3pvkzIxr3p3k1qr6SVbeMvibqvrytCfp7vu7e6m7l3bt2rWA2ABJtmCGmV/ALBZRyp5IcqCqrq2qK5LcluT4mjXHk3xo8hdMNyT5dXef7e5Pdvfe7t4/Oe8/uvuDC8gEMCszDBjCznkv0N3nq+quJI8k2ZHkS939bFXdOTl+X5ITSQ4nOZXk5SQfnvd5ARbBDANGUd1rPzoxvqWlpV5eXt7uGMAWqqonu3tpu3PMy/yC159Z55dv9AcAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAawkFJWVYeq6rmqOlVVx6Ycr6q6Z3L8maq6brJ/X1X9Z1WdrKpnq+qji8gDsBFmGDCCuUtZVe1Icm+Sm5McTHJ7VR1cs+zmJAcmt6NJvjDZfz7JP3b3nye5IcnfTzkXYNOYYcAoFvFK2fVJTnX3C939SpKHkxxZs+ZIkgd7xeNJrqyq3d19trufSpLu/m2Sk0n2LCATwKzMMGAIiyhle5K8uGr7dH5/KK27pqr2J3lXku8uIBPArMwwYAiLKGU1ZV9vZE1VvSnJV5N8rLt/M/VJqo5W1XJVLZ87d+6ywwKssekzzPwCZrGIUnY6yb5V23uTnJl1TVW9ISvD7KHu/trFnqS77+/upe5e2rVr1wJiAyTZghlmfgGzWEQpeyLJgaq6tqquSHJbkuNr1hxP8qHJXzDdkOTX3X22qirJvyU52d3/soAsABtlhgFD2DnvBbr7fFXdleSRJDuSfKm7n62qOyfH70tyIsnhJKeSvJzkw5PT353k75L8oKq+P9n3z919Yt5cALMww4BRVPfaj06Mb2lpqZeXl7c7BrCFqurJ7l7a7hzzMr/g9WfW+eUb/QEABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAFlLKqupQVT1XVaeq6tiU41VV90yOP1NV1816LsBmM8OAEcxdyqpqR5J7k9yc5GCS26vq4JplNyc5MLkdTfKFDZwLsGnMMGAUi3il7Pokp7r7he5+JcnDSY6sWXMkyYO94vEkV1bV7hnPBdhMZhgwhEWUsj1JXly1fXqyb5Y1s5wLsJnMMGAIiyhlNWVfz7hmlnNXLlB1tKqWq2r53LlzG4wIcFGbPsPML2AWiyhlp5PsW7W9N8mZGdfMcm6SpLvv7+6l7l7atWvX3KEBJjZ9hplfwCwWUcqeSHKgqq6tqiuS3Jbk+Jo1x5N8aPIXTDck+XV3n53xXIDNZIYBQ9g57wW6+3xV3ZXkkSQ7knypu5+tqjsnx+9LciLJ4SSnkryc5MOXOnfeTACzMsOAUVT31I9wDW1paamXl5e3Owawharqye5e2u4c8zK/4PVn1vnlG/0BAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAJQyAIABzFXKquqtVfVoVT0/uX/LRdYdqqrnqupUVR1btf+zVfXDqnqmqr5eVVfOkwdgI8wwYCTzvlJ2LMlj3X0gyWOT7Vepqh1J7k1yc5KDSW6vqoOTw48m+cvu/qskP0ryyTnzAGyEGQYMY95SdiTJA5PHDyT5wJQ11yc51d0vdPcrSR6enJfu/nZ3n5+sezzJ3jnzAGyEGQYMY95S9rbuPpskk/urp6zZk+TFVdunJ/vW+kiSb82ZB2AjzDBgGDvXW1BV30ny9imH7p7xOWrKvl7zHHcnOZ/koUvkOJrkaJJcc801Mz418Ho3wgwzv4BZrFvKuvu9FztWVb+oqt3dfbaqdid5acqy00n2rdrem+TMqmvckeSWJDd2d+ciuvv+JPcnydLS0kXXAaw2wgwzv4BZzPv25fEkd0we35HkG1PWPJHkQFVdW1VXJLltcl6q6lCSf0pya3e/PGcWgI0yw4BhzFvKPpXkpqp6PslNk+1U1Tuq6kSSTD4Ee1eSR5KcTPKV7n52cv6/Jnlzkker6vtVdd+ceQA2wgwDhrHu25eX0t2/SnLjlP1nkhxetX0iyYkp6/5snucHmIcZBozEN/oDAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAcxVyqrqrVX1aFU9P7l/y0XWHaqq56rqVFUdm3L841XVVXXVPHkANsIMA0Yy7ytlx5I81t0Hkjw22X6VqtqR5N4kNyc5mOT2qjq46vi+JDcl+Z85swBslBkGDGPeUnYkyQOTxw8k+cCUNdcnOdXdL3T3K0kenpx3weeSfCJJz5kFYKPMMGAY85ayt3X32SSZ3F89Zc2eJC+u2j492ZequjXJz7r76TlzAFwOMwwYxs71FlTVd5K8fcqhu2d8jpqyr6vqjZNrvG+mi1QdTXI0Sa655poZnxp4vRthhplfwCzWLWXd/d6LHauqX1TV7u4+W1W7k7w0ZdnpJPtWbe9NcibJO5Ncm+Tpqrqw/6mqur67fz4lx/1J7k+SpaUlbxMAMxlhhplfwCzmffvyeJI7Jo/vSPKNKWueSHKgqq6tqiuS3JbkeHf/oLuv7u793b0/K4PvummFDGCTmGHAMOYtZZ9KclNVPZ+Vvz76VJJU1Tuq6kSSdPf5JHcleSTJySRf6e5n53xegEUww4BhrPv25aV096+S3Dhl/5kkh1dtn0hyYp1r7Z8nC8BGmWHASHyjPwDAAJQyAIABKGUAAANQygAABqCUAQAMQCkDABiAUgYAMAClDABgAEoZAMAAlDIAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwACUMgCAAShlAAADUMoAAAaglAEADEApAwAYgFIGADAApQwAYABKGQDAAKq7tzvDhlXVuSQ/3e4cE1cl+eV2h7gMcm8tuef3J929a7tDzMv8Wgi5t5bc85tpfr0mS9lIqmq5u5e2O8dGyb215GZEr9Xfr9xbS+6t4+1LAIABKGUAAANQyuZ3/3YHuExyby25GdFr9fcr99aSe4v4TBkAwAC8UgYAMAClbB1V9daqerSqnp/cv+Ui6w5V1XNVdaqqjk05/vGq6qq6avNTz5+7qj5bVT+sqmeq6utVdeUm513v51dVdc/k+DNVdd2s546Yu6r2VdV/VtXJqnq2qj76Wsi96viOqvpeVX1z61JzOcwwM2wzcpthm6S73S5xS/KZJMcmj48l+fSUNTuS/DjJnya5IsnTSQ6uOr4vySNZ+W6iq14LuZO8L8nOyeNPTzt/gVkv+fObrDmc5FtJKskNSb4767mD5t6d5LrJ4zcn+dFrIfeq4/+Q5N+TfHMrMrvN9fs2w8ywzchthm3CzStl6zuS5IHJ4weSfGDKmuuTnOruF7r7lSQPT8674HNJPpFkKz/AN1fu7v52d5+frHs8yd5NzLrezy+T7Qd7xeNJrqyq3TOeO1zu7j7b3U8lSXf/NsnJJHtGz50kVbU3yfuTfHGL8jIfM8wMW3huM2xzKGXre1t3n02Syf3VU9bsSfLiqu3Tk32pqluT/Ky7n97soGvMlXuNj2Tl/zg2yyw5LrZm1n+GzTBP7v9TVfuTvCvJdxeecLp5c38+K/+B/t1mBWShzDAz7GLMsMHs3O4AI6iq7yR5+5RDd896iSn7uqreOLnG+y432yWfdJNyr3mOu5OcT/LQxtJtyLo5LrFmlnM3yzy5Vw5WvSnJV5N8rLt/s8Bsl3LZuavqliQvdfeTVfWehSfjsphhr2KGzc4MG4xSlqS733uxY1X1iwsv1U5e+nxpyrLTWfnMxQV7k5xJ8s4k1yZ5uqou7H+qqq7v7p8PnPvCNe5IckuSG3vyJvwmuWSOddZcMcO5m2We3KmqN2RlmD3U3V/bxJxrzZP7b5PcWlWHk/xhkj+qqi939wc3MS/rMMN+L/eFa5hhl2aGjTbDtvtDbaPfknw2r/6w6WemrNmZ5IWsDK8LHzr8iynrfpKt+5DsXLmTHEry30l2bUHWdX9+WXn/f/WHNv9rIz/7AXNXkgeTfH4b/p2+7Nxr1rwng31I1m3q79sM2/ysZtjW/jv9/3aGbXuA0W9J/jjJY0men9y/dbL/HUlOrFp3OCt/ffLjJHdf5FpbOdDmyp3kVFbej//+5HbfJuf9vRxJ7kxy5+RxJbl3cvwHSZY28rMfLXeSv87Ky+3PrPoZHx4995prDDfQ3Kb+rs0wM2zhuc2wzbn5Rn8AgAH460sAgAEoZQAAA1DKAAAGoJQBAAxAKQMAGIBSBgAwAKUMAGAAShkAwAD+F4YSboxTXC68AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "## QUESTION 4: Edit the code below to create a two-panel plot that shows r and v as a function of theta\n", "# Things to include:\n", "# - axis labels\n", "# - reasonable units for r and v, I suggest au and km/s\n", "# Check that your plot is reasonable! Are you getting values you expect for r and v?\n", "\n", "# This code creates two plots side-by-side, horizontally. You can do two vertical subplots using plt.subplots(2,1). \n", "# Then 'ax1.plot' is used instead of 'plt.plot' for the plot on the left, and ax2 is used for the plot on the right. \n", "# Figsize can be adjusted to make your plots bigger and space them out better. \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4))\n", "ax1.plot() # plot theta vs. r here in some units\n", "ax2.plot() # plot theta vs. v here in some units\n", "# Fill in between the '' to create labels below; see your previous python assignment\n", "ax1.set_xlabel('')\n", "ax2.set_xlabel('')\n", "ax1.set_ylabel('')\n", "ax2.set_ylabel('')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Caption your plot here: " ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 9)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m9\u001b[0m\n\u001b[0;31m ekin =\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "# QUESTION 5:\n", "# Next calculate the kinetic and potential energy of the system,\n", "# and the total energy. \n", "# We will plot these on the same plot to see how they vary. \n", "# You can do this using the equations you know for the kinetic and potential energy, \n", "# and by using the arrays for r and v that you calculated above you should then have\n", "# arrays for Ekin and Epot as a function of theta.\n", "\n", "ekin = \n", "epot = \n", "\n", "# Add together your ekin and epot arrays to calculate the total energy. Make sure they have\n", "# the same units! To test, you can add print statements to see what units you have in your \n", "# calculations, i.e.\n", "# Uncomment and edit this line:\n", "# print('The first value of ekin is {0}'.format(ekin[0].cgs))\n", "\n", "e_tot = \n", "\n", "# Lastly, calculate what Etot would be from the virial equation:\n", "e_tot_virial = \n", "\n", "# Check to see if your results agree! \n", "print('From the virial equation, E_tot = {0:.3}'.format(e_tot_virial.cgs))\n", "print('From the sum of the kinetic and potential energies, E_tot = {0:.3}'.format(e_tot.cgs[0]))" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "x and y must have same first dimension, but have shapes (10,) and (100,)", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# QUESTION 6: Plot your results!\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_rad\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mekin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'red'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'E_kin'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# edit this line to include units for ekin\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_rad\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mepot\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'blue'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'E_pot'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# edit this line to include units for epot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_rad\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0me_tot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcgs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'black'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'E_tot'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2809\u001b[0m return gca().plot(\n\u001b[1;32m 2810\u001b[0m *args, scalex=scalex, scaley=scaley, **({\"data\": data} if data\n\u001b[0;32m-> 2811\u001b[0;31m is not None else {}), **kwargs)\n\u001b[0m\u001b[1;32m 2812\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2813\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/__init__.py\u001b[0m in \u001b[0;36minner\u001b[0;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1808\u001b[0m \u001b[0;34m\"the Matplotlib list!)\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlabel_namer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1809\u001b[0m RuntimeWarning, stacklevel=2)\n\u001b[0;32m-> 1810\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1811\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1812\u001b[0m inner.__doc__ = _add_data_doc(inner.__doc__,\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(self, scalex, scaley, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1609\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmlines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLine2D\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_alias_map\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1610\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1611\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1612\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1613\u001b[0m \u001b[0mlines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_grab_next_args\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 391\u001b[0m \u001b[0mthis\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 392\u001b[0m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 393\u001b[0;31m \u001b[0;32myield\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_plot_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 394\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 395\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[0;34m(self, tup, kwargs)\u001b[0m\n\u001b[1;32m 368\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mindex_of\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtup\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 369\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 370\u001b[0;31m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_xy_from_xy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 371\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 372\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcommand\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'plot'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_xy_from_xy\u001b[0;34m(self, x, y)\u001b[0m\n\u001b[1;32m 229\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 230\u001b[0m raise ValueError(\"x and y must have same first dimension, but \"\n\u001b[0;32m--> 231\u001b[0;31m \"have shapes {} and {}\".format(x.shape, y.shape))\n\u001b[0m\u001b[1;32m 232\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 233\u001b[0m raise ValueError(\"x and y can be no greater than 2-D, but have \"\n", "\u001b[0;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (10,) and (100,)" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADUJJREFUeJzt3F9onfd9x/H3Z3YN6581oVFLZ6fUG05TXzQjUdMw1i1d2WpnF6bQi6SlYaFgwprSy4TB2ovcrBeDUpLUmGBCb+qLNbTuSBsGo80gSxcZUidOSNFcFmsuxGlLByksOPnu4pxNQpGtx+ccSY6+7xcI9JznJ+mrH/Lbj491nlQVkqTt73e2egBJ0uYw+JLUhMGXpCYMviQ1YfAlqQmDL0lNrBv8JMeSvJzkuYucT5JvJFlMcirJjbMfU5I0rSFX+I8ABy5x/iCwb/x2GPjm9GNJkmZt3eBX1RPAry6x5BDwrRp5CrgqyftnNaAkaTZ2zuBz7AbOrjheGj/2i9ULkxxm9K8A3vGOd9x0/fXXz+DLS1IfJ0+efKWq5ib52FkEP2s8tub9GqrqKHAUYH5+vhYWFmbw5SWpjyT/OenHzuK3dJaAa1cc7wHOzeDzSpJmaBbBPwHcOf5tnVuA31TVm57OkSRtrXWf0knybeBW4JokS8BXgbcBVNUR4DHgNmAR+C1w10YNK0ma3LrBr6o71jlfwBdnNpEkaUP4SltJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaGBT8JAeSvJhkMcl9a5x/d5LvJ/lpktNJ7pr9qJKkaawb/CQ7gAeBg8B+4I4k+1ct+yLwfFXdANwK/EOSXTOeVZI0hSFX+DcDi1V1pqpeA44Dh1atKeBdSQK8E/gVcGGmk0qSpjIk+LuBsyuOl8aPrfQA8GHgHPAs8OWqemP1J0pyOMlCkoXz589POLIkaRJDgp81HqtVx58CngF+H/gj4IEkv/emD6o6WlXzVTU/Nzd32cNKkiY3JPhLwLUrjvcwupJf6S7g0RpZBH4OXD+bESVJszAk+E8D+5LsHf9H7O3AiVVrXgI+CZDkfcCHgDOzHFSSNJ2d6y2oqgtJ7gEeB3YAx6rqdJK7x+ePAPcDjyR5ltFTQPdW1SsbOLck6TKtG3yAqnoMeGzVY0dWvH8O+MvZjiZJmiVfaStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJamJQ8JMcSPJiksUk911kza1JnklyOsmPZzumJGlaO9dbkGQH8CDwF8AS8HSSE1X1/Io1VwEPAQeq6qUk792ogSVJkxlyhX8zsFhVZ6rqNeA4cGjVms8Cj1bVSwBV9fJsx5QkTWtI8HcDZ1ccL40fW+k64OokP0pyMsmda32iJIeTLCRZOH/+/GQTS5ImMiT4WeOxWnW8E7gJ+CvgU8DfJbnuTR9UdbSq5qtqfm5u7rKHlSRNbt3n8Bld0V+74ngPcG6NNa9U1avAq0meAG4AfjaTKSVJUxtyhf80sC/J3iS7gNuBE6vWfA/4eJKdSd4OfAx4YbajSpKmse4VflVdSHIP8DiwAzhWVaeT3D0+f6SqXkjyQ+AU8AbwcFU9t5GDS5IuT6pWPx2/Oebn52thYWFLvrYkvVUlOVlV85N8rK+0laQmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqYlBwU9yIMmLSRaT3HeJdR9N8nqSz8xuREnSLKwb/CQ7gAeBg8B+4I4k+y+y7mvA47MeUpI0vSFX+DcDi1V1pqpeA44Dh9ZY9yXgO8DLM5xPkjQjQ4K/Gzi74nhp/Nj/S7Ib+DRw5FKfKMnhJAtJFs6fP3+5s0qSpjAk+FnjsVp1/HXg3qp6/VKfqKqOVtV8Vc3Pzc0NnVGSNAM7B6xZAq5dcbwHOLdqzTxwPAnANcBtSS5U1XdnMqUkaWpDgv80sC/JXuC/gNuBz65cUFV7/+/9JI8A/2TsJenKsm7wq+pCknsY/fbNDuBYVZ1Ocvf4/CWft5ckXRmGXOFTVY8Bj616bM3QV9VfTz+WJGnWfKWtJDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJamJQcFPciDJi0kWk9y3xvnPJTk1fnsyyQ2zH1WSNI11g59kB/AgcBDYD9yRZP+qZT8H/qyqPgLcDxyd9aCSpOkMucK/GVisqjNV9RpwHDi0ckFVPVlVvx4fPgXsme2YkqRpDQn+buDsiuOl8WMX8wXgB2udSHI4yUKShfPnzw+fUpI0tSHBzxqP1ZoLk08wCv69a52vqqNVNV9V83Nzc8OnlCRNbeeANUvAtSuO9wDnVi9K8hHgYeBgVf1yNuNJkmZlyBX+08C+JHuT7AJuB06sXJDkA8CjwOer6mezH1OSNK11r/Cr6kKSe4DHgR3Asao6neTu8fkjwFeA9wAPJQG4UFXzGze2JOlypWrNp+M33Pz8fC0sLGzJ15akt6okJye9oPaVtpLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDUxKPhJDiR5MclikvvWOJ8k3xifP5XkxtmPKkmaxrrBT7IDeBA4COwH7kiyf9Wyg8C+8dth4JsznlOSNKUhV/g3A4tVdaaqXgOOA4dWrTkEfKtGngKuSvL+Gc8qSZrCzgFrdgNnVxwvAR8bsGY38IuVi5IcZvQvAID/SfLcZU27fV0DvLLVQ1wh3Itl7sUy92LZhyb9wCHBzxqP1QRrqKqjwFGAJAtVNT/g62977sUy92KZe7HMvViWZGHSjx3ylM4ScO2K4z3AuQnWSJK20JDgPw3sS7I3yS7gduDEqjUngDvHv61zC/CbqvrF6k8kSdo66z6lU1UXktwDPA7sAI5V1ekkd4/PHwEeA24DFoHfAncN+NpHJ556+3EvlrkXy9yLZe7Fson3IlVveqpdkrQN+UpbSWrC4EtSExsefG/LsGzAXnxuvAenkjyZ5IatmHMzrLcXK9Z9NMnrST6zmfNtpiF7keTWJM8kOZ3kx5s942YZ8Gfk3Um+n+Sn470Y8v+FbzlJjiV5+WKvVZq4m1W1YW+M/pP3P4A/AHYBPwX2r1pzG/ADRr/Lfwvwk42caaveBu7FHwNXj98/2HkvVqz7F0a/FPCZrZ57C38urgKeBz4wPn7vVs+9hXvxt8DXxu/PAb8Cdm317BuwF38K3Ag8d5HzE3Vzo6/wvS3DsnX3oqqerKpfjw+fYvR6hu1oyM8FwJeA7wAvb+Zwm2zIXnwWeLSqXgKoqu26H0P2ooB3JQnwTkbBv7C5Y268qnqC0fd2MRN1c6ODf7FbLlzumu3gcr/PLzD6G3w7WncvkuwGPg0c2cS5tsKQn4vrgKuT/CjJySR3btp0m2vIXjwAfJjRCzufBb5cVW9sznhXlIm6OeTWCtOY2W0ZtoHB32eSTzAK/p9s6ERbZ8hefB24t6peH13MbVtD9mIncBPwSeB3gX9L8lRV/Wyjh9tkQ/biU8AzwJ8Dfwj8c5J/rar/3ujhrjATdXOjg+9tGZYN+j6TfAR4GDhYVb/cpNk225C9mAeOj2N/DXBbkgtV9d3NGXHTDP0z8kpVvQq8muQJ4AZguwV/yF7cBfx9jZ7IXkzyc+B64N83Z8QrxkTd3OindLwtw7J19yLJB4BHgc9vw6u3ldbdi6raW1UfrKoPAv8I/M02jD0M+zPyPeDjSXYmeTuju9W+sMlzboYhe/ESo3/pkOR9jO4ceWZTp7wyTNTNDb3Cr427LcNbzsC9+ArwHuCh8ZXthdqGdwgcuBctDNmLqnohyQ+BU8AbwMNVte1uLT7w5+J+4JEkzzJ6WuPeqtp2t01O8m3gVuCaJEvAV4G3wXTd9NYKktSEr7SVpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+Smvhf13t1pXPV3XUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# QUESTION 6: Plot your results! \n", "fig = plt.figure()\n", "# In the next three lines, change 'ekin' to ekin.to(some unit), etc. to make readable plots \n", "plt.plot(theta_rad,ekin,color='red',label='E_kin') \n", "plt.plot(theta_rad,epot,color='blue',label='E_pot') \n", "plt.plot(theta_rad,e_tot,color='black',label='E_tot')\n", "plt.xlabel('') # Add label\n", "plt.ylabel('') # Add label\n", "plt.legend(frameon=False)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using functions to explore how r and v vary with orbital parameters\n", "\n", "In python, once you start writing longer and more complicated programs, it is useful to create your own functions that you can then call in your code. Functions can be used for calculations that you might do many times, or to simplify a very long block of code. Remember you always have to run the function definition before you can run it in the code. This can be done by including your functions near the beginning of your notebook. Here, we will just define the function in the same block as running it. \n", "\n", "You create a function with a line starting 'def' , followed by your function name, brackets containing the parameters you need to give that function for it to do its work, followed by a colon. This tells python that everything underneath and indented (the indent is important!) that those lines are part of the function. \n", "\n", "We will just create a simple function here that returns two arrays, $r$ and $v$, but if you want to learn more one good tutorial is [here](http://introtopython.org/introducing_functions.html). " ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 5)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m5\u001b[0m\n\u001b[0;31m r =\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "# QUESTION 7: First create your function by copying and pasting your equations for r and v\n", "# below the 'def' statement in the space provided. Make sure you keep the indents! \n", "\n", "def calc_r_and_v(theta_rad, a, e, M):\n", " r = \n", " v = \n", " return r, v\n", "\n", "# We have told the function to return r and v arrays, given the inputs\n", "# The rest of the code is NOT indented because it is not part of the function! \n", "# You need to ensure that theta_rad is an array! There are tests you can do to ensure the \n", "# variable type is correct in a function like this, but we will leave that for another day.\n", "\n", "# Once you've defined the function, you can use it. Because our function returns TWO arrays,\n", "# we need to call it with two variables for them to go into:\n", "\n", "r_array, v_array = calc_r_and_v(theta_rad, a, e, M)\n", "\n", "# Print out a few numbers to check that you have arrays as outputs\n", "# and so you can see if your values match what you had earlier! \n", "r_array[0:5].to(u.au)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 13)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m13\u001b[0m\n\u001b[0;31m e1 =\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "# QUESTION 8: Calculate new r and v arrays for different eccentricity values, and plot together.\n", "\n", "# Now that we've defined a function, we can more quickly look at what happens when we change\n", "# some of the orbital parameters. \n", "\n", "# Use your function to plot r and v as a function of theta_rad as in Question 4, but add as\n", "# separate colours the values you get when changing e to a slightly larger value than you \n", "# used above for the Earth's orbit, and to a value close to, but less than 1. \n", "# Be sure to include a legend that states what e values you chose, and label your axes. \n", "\n", "# First, run your function to create two more arrays for r and v. We will keep all other\n", "# orbital parameters the same, but change e:\n", "e1 = \n", "r_e1, v_e1 = calc_r_and_v(theta_rad, a, e1, M)\n", "\n", "e2 = \n", "r_e2, v_e2 = calc_r_and_v(theta_rad, a, e2, M)\n", "\n", "# Create a plot\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4))\n", "\n", "# Plot the original results for Earth:\n", "# In these and following lines, be sure to change 'r' to r.to(some unit) and 'v' to v.to(some unit)\n", "# in order to make readable plots \n", "ax1.plot(theta_rad, r,color='black',label='') # plot theta vs. r here in some units\n", "ax2.plot(theta_rad, v,color='black',label='') # plot theta vs. v here in some units\n", "\n", "# Plot the results for your first eccentricity value:\n", "ax1.plot(theta_rad, r_e1,color='blue',label='') # plot theta vs. r here in some units\n", "ax2.plot(theta_rad, v_e1),color='blue',label='') # plot theta vs. v here in some units\n", "\n", "# Plot the results for your second eccentricity value:\n", "ax1.plot(theta_rad, r_e2,color='gray',label='') # plot theta vs. r here in some units\n", "ax2.plot(theta_rad, v_e2,color='gray',label='') # plot theta vs. v here in some units\n", "\n", "# Add appropriate labels below\n", "ax1.set_xlabel('')\n", "ax2.set_xlabel('')\n", "ax1.set_ylabel('') # include the units for r\n", "ax2.set_ylabel('') # include the units for v\n", "ax1.legend(frameon=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discussion\n", "\n", "Hopefully you can see that for the same semi-major axis value, and the same mass of the star and planet, changing the eccentricity of the orbit greatly changes values of $r$ at perihelion and aphelion. You should also see that when the planet is furthest from the sun, its speed in the orbit is slowest. If you don't see these things in your plot above, please re-check your calculations! \n", "\n", "### Question 9\n", "\n", "Do you expect the total energy in these systems has changed? Why or why not? Add your answer in the markdown text below (you could also use a scratch box to calculate this!):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Question 8 Answer:\n", "Your answer here! " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Extra box for stuff if you need it! " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }